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ABSTRACT 

A computer simulation model of the phosphate, phytoplank- 
ton and zooplankton dynamics in Monterey Bay was examined and 
modified. The model is driven by four forcing functions ex- 
pressed as annual cycles of upwelling velocity, incident 
solar radiation, mixed layer depth, and mixed layer tempera- 
ture. An alternate upwelling index was developed based on 
the local wind field. A revised radiation index is employed 
based on the generation of both advective fog and low stratus 
cloud cover common during upwelling on the California coast. 
Analysis of the model's response to sinking and advection of 
phytoplankton was examined. The importance of seasonal in- 
creases in predators was introduced as a controlling factor 
in the seasonal growth of zooplankton. The model is able to 
predict the seasonal trends of phosphate, phytoplankton, and 
zooplankton throughout the year. 
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I. INTRODUCTION 



A. PURPOSE 

In the interests of assessing the biomass of oceanic 
areas, various mathematical models have been designed to 
predict the dynamic response of ecosystems to change in the 
physical and chemical environment. Such models necessarily 
rely on accurate characterization of the forcing functions 
and accurate representation of the system with equations 
which are consistent with the specific region under study. 

Recent attention paid to modeling systems in upwelling 
regions, Coastal Upwelling Ecosystems Analysis (CUEA),has 
brought about developments in both these areas. Adequate 
data of a physical and chemical nature is becoming available 
and significant progress in the refinement of governing equa 
tions is evident (Walsh, 1973). 

The research presented here was aimed at creating a simu 
lation model to describe the seasonal plankton dynamics in 
Monterey 3ay, California, as known from the best available 
data. An existing simulation model (Pearson, 1975) was 
evaluated in an effort to make refinements and to judge the 
applicability of the time simulation technique to the 
Monterey region, characterized by seasonal upwelling. 

3. MODELING PHILOSOPHY IN THE ECOLOGICAL CONTEXT 

The merit of mathematically modeling systems of a purely 
physical nature in which parameters, initial conditions, 
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boundary conditions and the variable relationships are readi- 
ly specified is unique. The value of current models of a 
biological nature requires some comment. 

According to Odum (1971) there are a number of reasons 
behind constructing any mathematical model. Prediction of 
future states of the variables involved is an end in itself. 
Use of the prediction accuracy to focus data acquisition 
efforts or direct attention to deficiencies in concepts is 
likewise valid. If the model is real in the sense that it 
attempts to predict changes in the state variables based on 
a framework drawn from research in the real world, then 
failures may be useful in delineating weaknesses in the 
governing equations. 

Construction of a "real" model is at best difficult. A 
perfect model of an ecosystem assumes that "all 'states and 
rates of change of the variables are known at all times" 
(Walsh, 1973). Given the vastness of complex interactions 
in even the simplest of biological systems , it may be stated 
that "no perfect representation of the real world exists 
except the real world itself" (Walsh, 1973). 

It follows then that certain assumptions, logical state- 
ments and approximations of real dynamics are necessary to 
allow the model to operate. The resulting model may bear 
little correlation to the actual physical, chemical and 
biological situation at hand. Again, the argument for con- 
structing the model must be considered. For purposes of 
prediction, a model which makes a logical statement relating 
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two state variables in the attempt to specify a relationship 
which may be unknown in part or in whole has merit. The 
fact that the logic of the model is a crude simplification 
is incidental if the prediction capability is proven. 

All models can be characterized on the basis of realism, 
precision, and generality (Odum, 1971). It has been suggested 
that these three qualities are mutually exclusive in ecologi- 
cal models (Patten, 1971). A fourth category, simplicity, 
might be considered. Exacting detail, while lending to 
realism and precision often forces the model to be specific. 
Only when the model becomes an isomorph, i.e., a one-to-one 
correspondence to the real world, will generality be restored 
(Walsh, 1973). 

The problems in creating a realistic model have been dis- 
cussed briefly. The value of this type of model lies in re- 
search guidance. Where a specific question is asked of an 
ecological model, precision in prediction capability may be 
enhanced by sacrificing the other qualities. Such an 
approach might be taken in a fisheries model where a precise 
output for a limited region is desired (Odum, 1971). Simplic- 
ity and generality at the expense of both precision and 
realism may enhance understanding of broad ecological con- 
cepts . 

The study presented herein was conducted with the idea 
of creating a model having good predictive capabilities for 
conditions occuring over a long span of time in a limited 
oceanic region. 
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C. DESCRIPTION OF THE MODEL 



The original ecosystem model was developed by Pearson 
(1975) as part of continuing research into the correlation 
of zooplankton biomass with the chemical and acoustic pro- 
perties of the ocean (Traganza, 1976). It is the intent of 
this research to identify areas of weakness in the model 
and make the appropriate changes to the equations and input 
parameters. These changes include a refinement of the up- 
welling index and radiation index as well as inclusion of 
advection, sinking and predation terms in the governing 
equations . 

The modeling technique used is a time dependent simula- 
tion solved digitally on an I3M 360 computer using the IBM 
Continuous Simulation Modeling Program (CSMP-360). The 
dynamic equations are written in non-linear differential 
form and driven by four exogenous forcing functions expressed 
as annual cycles of upwelling velocity, mixed layer depth, 
mixed layer temperature, and incident solar radiation. The 
output consists of annual cycles of phosphate concentration 

(.micro gram- atoms of phosphorous per liter, ug-atP/Jl), phyto- 

2 

plankton biomass (grams of carbon per square meter, gC/m ) 

2 

and standing stock of herbivorous zooplankton (g C/m ) in 
the mixed layer. 

The Pearson version (Figure 1) of the model is defined 
by the following basic equations: 

1) = NUT + REGEN - UPTAK (1) 
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LIGHT 




Figure 1. Pictoral representation of Monterey Upwelling 
Ecosystem. 
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where : 



gg— = time rate of change of phosphate concentra- 
tion (yg-at P/1) in the mixed layer 

NUT = input of phosphate from upwelling and mixing 
from below the mixed layer (yg-at P/1) 



REGEN = phosphate recycled as biological excretory 
products within the mixed layer (yg-at P/1) 

UPTAK = phosphate depleted by phytoplankton utiliza- 
tion (yg-at P/1) 



2 ) 



dX2 

dt 



where : 



dX2 

dt 



3) 



PROD 

RESP2 

GRAZ 

dX3 

dt 



PROD - RESP2 - GRAZ (2) 

time rate of change of phytoplankton biomass 
( gC/m~ ) 

2 

photosynthetic production (gC/m ) 

2 

respiration losses (gC/m ) 

2 

grazing losses (gC/m ) 

GRAZ - RESP3 - VOID - LOSS (3) 



GRAZ = 

RESP3 = 
VOID = 
LOSS = 



input 9 due to ingestion of phytoplankton 
(gC/ni ) 

2 

respiration losses (gC/m ) 

2 

excretory losses (gC/m ) 

2 

stock depletion by mortalxty (gC/m ) 

(Pearson, 1975) 
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II. BACKGROUND THEORY 



A. UPWELLING THEORY 
1 . Ekman Dynamics 

Seasonal upwelling of nutrient-rich deep water into 
the productive mixed layer is an important aspect of the 
-Monterey Bay ecosystem. The principal source of phosphate 
in the mixed layer is vertical advection which is produced 
in response to wind-induced Ekman circulation. Theoretical 
calculations of the offshore component of horizontal current 
are based on Ekman pure drift theory in an infinitely deep, 
homogeneous ocean (Neumann and Pierson, 1966). The compo- 
nents of the horizontal current take the form: 



-(tt/D)Z 

V. e 
o 


COS [45- ( 7T / D ) Z ] and 


(4a) 


,, - ( tt/D) Z 

V o e 


SIN [45-(tt/D)Z] 


(4b) 


(A/pa) SIN 4) ) 


is 

and 


(5a) 


t /Aa 




(5b) 



o y 

where : 

U = velocity component in the x direction 

V = velocity component in the y direction 
(parallel to wind) 

V = surface velocity component (45° to right of 
wind, Northern Hemisphere) 

Z = depth 

D = Ekman depth of frictional resistance 
A = coefficient of eddy viscosity 
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Ty = surface wind stress in y direction (dynes) 
a = fp/A 

<p = latitude 

. 3 

p = density of water (g/cm ) 

f = Coriolis parameter = 2 to SIN <j> 

The coefficient of eddy viscosity, A, was assumed 
constant since detailed information on the small scale mo- 
tions of the surface layer was not available. This assump- 
tion is generally made in mass transport studies (Neumann, 
1968) . 

The mass transport due to wind driven current is 
found by integrating the equations of velocity (U and V) 
over the depth of the water column, as shown: 

S = p UdZ = t /f (6a) 

x y 

and S = p VdZ = 0 (6b) 

y 

where: 

= mass transport in the x direction 
Sy = mass transport in the y direction 

From the equations for mass transport, it is observed that 
the net transport is 90 degrees to the right of the wind in 
the northern hemisphere and that it is directly proportional 
to the wind stress acting on the sea surface, parallel to 
the coast. 
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2. Wind Stress on the Sea Surface 



The coupling of wind energy to the water at the air- 
sea interface is defined by the conventionally-accepted form 
of the wind stress equation: 

T = p'C W 2 (Gordon, 1972) (7) 

y £a La 

2 

where: = wind stress parallel to the coast (dynes/cm ) 

p' = density of air 

C = drag coefficient at height z (non-dimensional) 

W = wind speed at height z (cm/sec) 

La 

As in other momentum exchange studies, the dynamic 
coefficients appear to present the singular most difficult 
problem. Wilson (1960) states that the drag coefficient for 
air flow over a water surface is dependent on the wind speed 
(W) , the height of measurement ( z ) , a surface roughness para- 
meter (z^), and atmospheric stability terms. There may be 
additional dependence on oceanic parameters of depth, fetch 
distance and wave height. An average value of the drag co- 
efficient for high wind speeds (i.e. greater than 10 m/sec) 

- 3 . 

of 2.37 X 10 is arrived at by examination of research com- 
pleted through 1959 (Wilson, 1960). 

Work by Deacon (1962) established an empirical rela- 
tionship for winds up to 13 meters/second which gives a 
linear dependence of C z (at z = 10m) to the wind speed as 
follows : 

C 1Q = (1.10 + 0.04 W 1Q ) X 10 -3 for W in m/sec (8) 

(Roll, 1965) 
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Deacon’s research was carried out using data from ship obser- 
vations and coastal regions . 

3 . Spatial Extent 

The offshore movement of water in a region bounded 
by a coastline causes replacement water to be upwelled from 
below the layer affected by the surface wind stress. Accord- 
ingly, it is necessary to specify the spatial dimensions of 
the region in order to apply principles of continuity in de- 
termining the vertical current speed. The significant para- 
meters involved in upwelling are: (1) the offshore horizontal 
dimension of the upwelling region, which is most directly 
related to the width of the wind field; (2) the offshore 
component of the Ekman current (U); and (3) the depth of the 
Ekman layer (D). 

Estimates of the maximum depth of the source of up- 
welled water in a coastal region are approximately 100 to 
200 meters (Neuman and Pierson, 1966) based on the slope of 
the isotherms. The mass transport offshore due to wind 
drift current has been previously detailed as occurring in 
a surface layer of depth, D, corresponding to the Ekman depth 
of frictional resistance; regardless of the source depth of 
the upwelled water, it transgresses the "boundary" at depth, 
D, before being carried offshore. The vertical dimension of 
the region is then readily specified. It is assumed that 
the depth of frictional resistance (D) and the depth of the 
mixed layer (Z) are approximately coincident. The differ- 
ences between (D) and (Z) are important when phytoplankton 
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are advected out of the Ekman layer but exist throughout 
the mixed layer. 

The horizontal extent of the upwelling region is 
considerably more nebulous than the depth of the mixed layer. 
Classical estimates of this dimension limit the zone to 
probably no more than 100 kilometers in offshore width. It 
may be expected that the width of the region is dependent 
on several factors; among them, the width of the wind field, 
variations and non-linearities in the energy exchange pro- 
cesses due to surface roughness characteristics or stability 
changes in both the atmosphere and the water column, as well 
as spatial and temporal oscillations of the significant wind 
vector. In addition, local coastal and bottom topography 
may figure extensively in the problem as will be discussed 
later . 

Sverdrup (1938) observed a relatively well defined 
offshore boundary to a coastal upwelling zone, coincident 
witha downwind current which was marked by an intense verti- 
cal gradient of velocity (Sverdrup et al., 194-2). He further 
observed an offshore migration of the current band at a rate 
somewhat less than the speed of the offshore surface current 
within the upwelling region. The latter fact implied the 
possibility of cellular circulation patterns in the near sur- 
face upwelling zone. 

Hidaka (1954) proposed a steady state theory of up- 
welling in which he defined a horizontal frictional distance, 
analogous to the depth of frictional resistance 
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appearing in Zkman 1 s work. Applying reasonable values to 
Hidaka’s expression confines the total cellular circulation 
region to a width of about 15 kilometers. The width of the 
region is defined by: 

°H = ^2A h /2w SIN where: A R = 10 7 (gcm' 1 sec -1 ) (9) 

(Smith, 1968) 

Downward vertical currents occurring in the seaward half of 
the cell limit the upwelling to a width of 7.5 kilometers. 
The theoretical width is not consistent with observations of 
upwelling at greater distances from the coast (Barnes, 1969) 
Yoshida (1967) developed a quasi-steady state model 
applicable to eastern boundary current regions (Smith, 19 6 S' 1 
An expression was derived for the horizontal dimension of 
the coastal upwelling region which is given by: 

L, = - 1 1 — (Hy/Dy)^ for latitudes greater than 

22 degrees (10) 

where: L = horizontal width of the coastal boundary 

region (m) 

_2 

g = gravitational acceleration (m-sec ) 

H = thickness of the upper layer (m) 

D = H + thickness of the lower layer (m) 

f = Coriolis parameter (sec -1 ) 

_ 3 

p = density of water (g-cm ) 

Ap = density difference between deep and surface 
layers (g-cm~ 3 ) 

y = internal friction coefficient (sec -1 ) 
y = 2 X 10 -7 sec -1 (Smith, 1968) 
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3. CURRENTLY AVAILABLE UPWELLING INDEX 



Bakun (1976) has calculated coastal upwelling indices 
for the coast of North America using Fleet Numerical Weather 
Central monthly mean pressure fields to compute the geos- 
trophic wind. Analysis is done on a 63 x 63 point grid of 

approximately 200 nautical mile mesh length (Bakun, 1973). 

/ 

Figure 2 shows Bakun's results for 1974 at 36 degrees North, 
122 degrees West, approximatley 50 miles south of Monterey 
Bay. In reviewing several aspects of hydrographic surveys 
taken in Monterey Bay, it was suspected that this region is 
characterized by local upwelling patterns not appearing in 
the index computed by Bakun. Further discussion follows in 
the methods section. 

C. EFFECTS OF UPWELLING ON INCIDENT RADIATION 

Cold water upwelled from depth during the spring and 
summer months brings about both advection fog and stratus 
cloud formation. The fog occurs as warm surface air is 
cooled by the cold seawater to its saturation point causing 
the condensation of the contained water vapor. Stratus 
clouds occur below the base of a quasi-permanent atmospheric 
temperature inversion (Tont, 1975). 

The net effect of both fog and stratus layers is to re- 
duce the amount of solar radiation reaching the sea surface 
during the upwelling period. Tont (1975) obtained mean 
values of solar radiation at the surface over the period 
1950-1973 at San Diego. A correlation study of the annual 
upwelling index and the amount of sunlight at the surface 
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Figure 2. Upwelling index for 1974 at 36 degrees North, 
122 degrees West, computed by Bakun. 



showed a maximum of about 40 per cent reduction in available 
sunlight during the peak upwelling months. The result of 
Tont's research is shown in Figure 3. 

D. SINKING OF ALGAL CELLS 

Review of sedimentation rates of algal cells indicate 
justification for including a sinking term in the phytoplank- 
ton equation. The Pearson (1975) model assumes uniform algal 
concentation within the mixed layer. It is reasonable that 
changes in the concentration of phytoplankton caused by ver- 
tical movement out of the layer will occur. 

Parsons and Takahashi (1973) state that the theoretical 
sinking rate of phytoplankton can be determined from: 

V = (fi-^fi.) (11) 

9 n*<P r 

where : 

r = radius of the cell 

g = gravitational acceleration 

p' = density of the organism 

p => density of the medium (sea water) 

n = viscosity of the medium 

<P = form resistance coefficient (accounts for non- 
spherical shapes) 

Measured rates of sinking for live algal cells vary be- 
tween zero and 30 meters per day according to Parsons and 
Takahashi (1973). These rates represent a broad range of 
species and sizes. 
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Figure 3. Monthly means of upwelling index (m /se 
and percent of possible sunshine (after 
Tont, 1975). 



ADVECTION 



Cj . 

It is known that significant current velocities occur 
within the mixed layer from Ekman dynamics. It might be ex- 
pected then that the concentration of phosphate, phytoplank- 
ton and zooplankton will be affected by advection. O’Brien 
and Wrobleski (1972) have shown by scale analysis that advec- 
tion and biological productivity of the Peru upwelling eco- 
system are equally important. 

The advective terms of the material derivative of a 
property (C) are written: 

Advective change = U + V 1^- + W -|S- (12) 



where: U,V,W = velocity components in x, y, and z directions 

respectively 

3C 3C 3C 

ttt, -ry = gradients of the concentration of property 

(C) in the x, y, and z directions 



A review of California Cooperative Fisheries Investiga- 
tion (CALCOFI) data for 1974 (CALCOFI, 1974) shows that the 
phosphate concentration varies significantly in the vertical 
direction but the horizontal variation is not as significant. 

It is assumed for the purposes of this model that the ad- 
vection term for phosphate reduces to the vertical component: 



Advective change = W 



9X1 



since : 



9Z 

9X1 9X1 



(13) 



= 0 



9X 9Y 

Pearson (1975) has included this term in the computer model 

DX1 - XI 



as : 



UPWELL = W 



EKL 



(14) 
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where : 



W = vertical speed (m/day) 

EKL = depth of the Ekman layer (m) , EKL = D 

DX1 = concentration of phosphate below the mixed 
layer (yg-at P/1) 

XI = concentration of phosphate in the mixed layer 
(constant) (yg-at P/1) 



DX1 - XI 
EKL 



vertical gradient of phosphate (yg-at P/m 1) 



The advection of zooplankton in the revised model is 
discounted because it is assumed that they have developed 
mechanisms which permit them to remain in a particular 
oceanic region by swimming or riding cellular circulation 
currents . 

Phytoplankton, however, are affected by advection. 

CALCOFI data (CALCOFI, 1974) shows that the significant 
gradients of phytoplankton occur in the vertical and off- 
shore directions; therefore, the advection expression for 
phytoplankton is reduced to: 

Advective change = U + W (15) 

The horizontal velocity is related to the vertical veloc- 
ity by continuity (when it is assumed that = parallel 
to the coast) such that: 

U = - WL/D (16) 
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where : 



U = horizontal velocity component offshore (m/day) 
W = upward vertical velocity component (m/day) 

L = width of the region (m) 

D = Ekman depth (m) 

The phytoplankton advection term can be written: 



Advective change = W [- ^ 



(17) 



The gradients of phytoplankton can be approximated by: 



3X2 LIM AX2 ^ 3X2 LIM AX2 

3X AX+0 AX anQ 3Z ‘ AZ+0 AZ 



(18a6b) 



If it is assumed that the concentrations of phytoplankton 
below the mixed layer and outside the modeled area are much 
less than the concentration in the modeled region, then the 
gradient can be approximated by: 



X2-0 

AX 



horizontal gradient; 



(19) 



X2-0 

AZ 



vertical gradient; 



( 20 ) 



where: X2 = phytoplankton concentration in the 

mixed layer 

AX = horizontal distance at which the concen- 
tration of phytoplankton becomes negli- 
gible with respect to the modeled area 

AZ = vertical distance at which phytoplankton 
concentration becomes negligibly small 

The phytoplankton advection term can then be written as: 

Advective change to T , 

Phytoplankton in the = W(X2)[- + -gr) (21) 

modeled area 
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The sign convention is established as positive (+) into the 
model area and negative (-) out. 

F. PREDATION PRESSURE 

Under constant environmental conditions, the dynamics of 
zooplankton growth are a function of the food supply and 
predation pressure. Natural processes of population control 
insure that a given species will not be eliminated by low 
food supply and high predation pressure or that unstable 
growth will not result from the inverse situation. Attempts 
to translate the stabilizing factors into mathematical terms 
often result in simplified general relationships that do not 
exhibit the flexibilities inherent in the real system. Not- 
withstanding the limitations, a predation term is included 
in the simulation model. 

The long term survivability of both predator and prey is 
keyed directly to maintenance of a balance in energy expendi- 
tures and gains. To date, laboratory experiments to dupli- 
cate this balance have not been entirely successful. Prey 
stocks in small laboratory environments have been artificial- 
ly supported by providing refuge where predator access to 
the prey was denied and by introducing additional prey to 
the system to replenish the supply (Patten, 1971). Neverthe- 
less, population control by predation is prey-density depen- 
dent as a first approximation. For a system of one prey 
species and a single predator species, predation pressure 
may be represented by: 

P = k^ D for 0<D<d^ and 

P = k 2 for D>d 1 
29 



( 22 ) 



where: P = predation pressure or food intake per unit 

predator 

D = prey density 

d^= satiation threshold of prey density 
k^ and k 2 = constants 

Above prey density values of d^, predator satiation will oc- 
cur and pressure will level off at some constant high value. 
The function is shown on an arbitrary scale in Figure 4 
according to Patten (1971). This simple relationship is 
not entirely satisfactory for the predator since low food 
supplies would result in starvation. Under such circum- 
stances, and given a second prey species of perhaps less 
palatability , the predator would switch food intake to the 
more abundant species of prey. There are several obvious 
implications to the switching phenomenon: 

1. predation on the most abundant species assures 
stable growth of that species , 

2. switching from a declining species serves to prevent 
elimination of that species (Patten, 1971), and 

3. elimination of the predator by overgrazing one avail- 
able food source is avoided. 

The model under consideration in this thesis is a single 
species approximation in which predation pressure is incor- 
porated with other zooplankton mortality factors, i.e. 
morphological death, in the "LOSS" term. This term is given 
by Pearson (1975) as: 

LOSS = L (X3 ) (23) 
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PREDATION PRESSURE 




Figure 4. Predation pressure as a function of prey 
density (zooplankton biomass) (after 
Patten , 1971) . 
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where : 



LOSS = fraction of the zooplankton population re- 
moved from the system in a day’s time 
(g C/m 2 day) 

L = constant rate of loss or percent loss per 
day (day“l) 

2 

X3 = zooplankton "standing crop" (g C/m ) 

Harriston concluded in 1960 that "herbivores are preda- 
tor limited" (Patten, 1971). 

Conclusions by Patten (1971) show that herbivores are 
both food and predator limited. Therefore, in the modified 
model, the LOSS term is reserved for variable predator limi- 
tation of herbivorous zooplankton and natural death of 
zooplankton herbivores is assumed to be negligible. 
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III. METHODS 



The existing model appears weak in several areas. Orig- 
inally, the forcing functions of upwelling and solar radia- 
tion which are critical to the simulation were obtained from 
average cycles for the California coast (upwelling) and the 
latitude band between 30 and 40 degrees north. Phytoplankton 
sinking and advection terms were lacking and the predation 
pressure term in the zooplankton equation needed improvement. 

A. UPWELLING PROGRAM 

A computer program was written to calculate the annual 
upwelling index for Monterey Bay. Wind data was obtained 
from the Monterey Peninsula Airport observations during 1974. 
Hourly observations were examined to determine the signifi- 
cant mean daily wind vector for each day and corrections 
were applied to speed and direction when indicated by topo- 
graphic obstructions. Fortunately, no corrections were 
needed to northerly winds which are the driving mechanism 
for upwelling. 

An x-y coordinate system was established such that the 

2 

surface wind stress (t = p'C W ) could be directly calcu- 
lated from the wind component parallel to the coast. The 
Ekman mass transport (S ) , surface current velocity (V Q ) 
and the average velocity within the Ekman layer (U) were 
calculated by the following relations: 

V = t /Aa (5b) 

o y 
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where : V = 

o 

and 

U = 

where: U = 

and : 

S = 
x 

where S = 

x 

The vertical current representing upwelling is determined 
per unit y (along coast) from the continuity equation, 

A 

(V • v) = 0, such that 

0(D) = W(L) (25) 

where: U _= average horizontal velocity in the Ekman 

layer (m/day) 

W = upwelling speed (m/day) 

D = Ekman depth (m) 

L = horizontal width of the upwelling region (m) 
Zero current is assumed parallel to the coast. 

The horizontal width of the upwelling region, (L) was 
determined from Yoshida's (1967) equation. Applying typical 
values to this equation yields a horizontal width of 3 X 

4 

10 m. 

Principles of continuity were used to arrive at the mean 
daily upwelling current values and these figures were aver- 
aged over a 30-day period. 



surface current vector 



*/ 



UdZ 



( 24 ) 



average velocity in x direction in the layer 



x / f 

y 



(6a) 



mass transport in x direction (offshore) 
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B. RADIATION INDEX 



A revised incident solar radiation index which takes in- 
to account the reduction of sunlight by local fog and stratus 
cloud cover was developed. Correlation of the measurements 
of surface irradiance and the strength of upwelling (as in- 
dicated by the upwelling index) by Tont (1975) is repre- 
sented by Figure 5. This figure is derived from the results 
of Tont ' s study as shown in Figure 3 by plotting the percent 
of possible sunshine (Y-axis) as a function of the upwelling 
index (X-axis) which has been scaled to a maximum value of 
one. The resulting curves (A) and (B) show the conditions 
occurring after, and prior to, the upwelling maximum, respec- 
tively. 

The difference exhibited between the curves may be due 
to changes in the character of the air mass brought about by 
seasonal migrations of the quasi-permanent thermal low over 
Nevada and the North Pacific sub-tropical high located west 
of the coast. 

Figure 5 was used to calculate a revised incident radia- 
tion index by entering with the newly-developed upwelling 
index values and by applying the corresponding percent to 
the theoretical mean monthly radiation for clear sky condi- 
tions (possible sunlight). Figure 6 shows the possible sun- 
light throughout the year (from Tont, 1975) as Q , and the 
theoretical radiation at the surface (Qk) after the effects 
of fog and clouds are considered. 
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PERCENT SUNLIGHT 



100 



90 - 




Figure 5. 



Percent of possible sunlight as a function 
of upwelling index. 
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Figure 6. Seasonal variation of incident solar radiation, clear sky 
conditions (Q ) and reduced by fog and clouds (Q-). 



The incident solar radiation is an essential part of 
the ecosystem simulation model in that it is used in the 
equations governing the production of organic carbon by the 
photosynthetic activity of phytoplankton. These equations 
are discussed in detail by Pearson (1975). 

C. ALGAL SINKING TERM 

The phytoplankton sinking term which was added to the 

Pearson (1975) model is written in the form used by Riley 

(1965) in the units of the quantity of phytoplankton trans- 

2 

ferred per day (g C/m -day). The amount of phytoplankton 
that sinks out of the mixed layer in a day is given by: 

SINK = ( V/Z ) X2 (26) 

where: SINK = flux of phytoplankton out of mixed layer 

(gC/m 2 -day ) 

V = ^sinking rate (m/day) 

Z = mixed layer depth (m) 

2 

X2 = phytoplankton biomass in mixed layer (gC/m ) 

V/Z = fraction of phytoplankton which sinks out of 
mixed layer per day 

Vertical circulation with the water column becomes 
significant as the upwelling current speed and sinking rates 
of phytoplankton approach the same order of magnitude. Since 
the phytoplankton are non-mobile and are generally of the 
same density as the water or have developed shapes which 
increase their sinking resistance, they will be carried along 
with vertical currents in the water column. The vertical 
circulation during upwelling opposes sinking but may 
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accelerate downward transfer during brief periods of surface 
convergence and associated downwelling. 

The complete phytoplankton sinking term is: 

SINK = ((V - W)/Z)X2 (27) 



where: W = upwelling speed (m/day) 

An average value of one meter per day was used 'for the sink- 
ing rate (V) in the computer simulation model. This value 
was determined from estimates presented by Lehman et al . 
(1975) and Bannister (1974). The SINK term acts to decrease 
the phytoplankton concentration in the modeled region and is 
subtracted in the phytoplankton equation (2) as shown: 



d X2 
dt 



PROD 



RESP2 



GRAZ 



SINK 



(28) 



D. PHYTOPLANKTON ADVECTION TERM 

An equation describing the horizontal advection of phyto 
plankton is included in the revised model. This term is in 
the units of flux and is written as : 



ADVEC2 = W (X2)K (29) 



where: • ADVEC2 = the change of concentration of phytoplank 

ton over time (in the mixed layer) due to 
advection (g C/m2-day) 

W = upwelling speed (m/day) 

X2 = phytoplankton concentration in the mixed 
layer (g C/m^) 

K = advection coefficient (m ^) 
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The simplified advection term used in the model is de- 
rived from the equation (21) shown in the background theory 
section, which is: 

Advective change = W (X2) [- 9-r^- + tt] (21) 

The Ekman current offshore which results in an upwelling 
current (W) is effective over the depth of frictional resist- 
ance (D) , but the phytoplankton in the model (which are 
uniformly distributed in the mixed layer) will be advected 
at the average velocity over the depth of the mixed layer as 
determined by the fraction of the mixed layer that is coinci- 
dent with the Ekman layer. It is known that the mixed layer 
varies seasonally from about 10 to 100 meters in depth and 
this causes an annual (seasonal) variation in the average 
horizontal current in the phytoplankton advection term. The 
average current of the mixed layer can be determined from: 

z 

U = =r J' U„dZ where: U^. = a function of e 

o 

(30) 

A shallow mixed layer experiences higher average veloci- 
ties than a deep layer as seen from the integral expression 
(eq. 30 ) . 

A coefficient to describe the seasonal average of the 
vertically averaged horizontal current in the mixed layer is 
included in the advection equation as shown: 

Advective change = W (X2) [- + (31) 
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The coefficient, K , represents the seasonal average of the 
fraction of the mixed layer which is coincident with the 
Ekman layer and is affected by the Ekman current. 

The constants in the advection expression (within the 
brackets of eq. 31.) are lumped into the advection coefficient, 

K, for convenience in equation (29). The range of K varies 

-1 -2 -1 

rrom about 10 to 10 m depending on the estimated values 
of the gradient distance, (AX) and (AZ) and the velocity co- 
efficient, K . A value of 0.1 m ^ was used in this model 
m 

for K. 

E. PREDATION LOSS 

The predation pressure function hypothesized for Monterey 
Bay is based on the assumption that the system supports a 
resident population of herbivorous zooplankton predators 
throughout the year. The pressure exerted on the herbivore 
prey by this maintenance level (N) of predators follows the 
curve in Figure 7 below the prey density threshold, d. The 
pressure is sufficiently low to allow growth of the zooplank- 
ton stock, but high enough to stabilize growth. A second 
pressure is imposed on the zooplankton (see Figure 7) after 
their density reaches the critical level, d. The simulation 
approximates conditions which may be brought about as tran- 
sient predators move into the area presumably in response to 
increased food availability. In summary, predation pressure 
may be represented by a two part function with a pressure- 
density curve, a, due to resident predator species and curve 
b due to the transient predator population. 
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PREDATION PRESSURE 



d 




l l l 

PREY DENSITY 



Figure 7. Predation pressure on zooplankton by 
transient predators. 
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Predator biomass may be expressed as shown in Figure 8, 
as a simple step increase or some other function of prey 
density. Consideration of the step type function is impor- 
tant because it allows the predator population (and conse- 
quent pressure) to maximize during those seasons when the 
zooplankton "standing crop" will support additional numbers 
of predators; the high predator pressure will deplete the 
food stock to a lower level than would be reached by the 
resident predator pressure alone. Measurements of zooplank- 
ton biomass for 1974 (Traganza, 1975) suggest a rapid decline 
in "standing crop" following an early summer maximum. One 
might suspect that a transient predator population is forced 
out of the region once the food sources are depleted. This 
is done in the model by decreasing the predator population 
when the zooplankton biomass declines to a specified level. 
The predator levels were related directly to zooplankton 
"standing crop" by triggering an increase in predator popula- 
tion at a threshold of zooplankton biomass during periods 
when this prey population was on the upswing. The predators 
were similarly reduced at a second threshold during the de- 
dining phase of zooplankton growth. A set of four condi- 
tional statements is used in the simulation to describe the 
predator- prey relationship. These statements specify the 
following : 



dX3 

dt 


0 


AND 


X3 


< 1.0 


THEN : 


FISH = 1.0 


dX3 > 
dt 


0 


AND 


X3 


>1.0 


THEN: 


FISH = 5.0 
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PREDATOR POPULATION 



d 




Figure 3. Predator population. 
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0 AND X3 > 0.2 THEN: FISH = 5.0 

0 AND X3 < 0.2 THEN: FISH =1.0 

The modified predation pressure term in the zooplankton 
equation of the model is: 

LOSS = (L) (FISH) (X3) (32) 



IF 



dX3 
dt — 



IF 



dX3 

dt — 



where : 

LOSS = the amount of herbivorous zooplankton biomass 
lost per day as a result of predation 
(g C/m^ day) 

L = loss rate or percent loss per FISH per day 
(day"l) 

FISH = number of predators (non-dimensional) 

2 

X3 = zooplankton "standing crop" (g C/m ) 
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IV. RESULTS 



A. UPWELLING CALCULATIONS 

The mesoscale wind data used to calculate upwelling for 
this simulation is shown by Figure 9. The plot depicts the 
frequency of occurrence (f = N/EN) for 45-degree segments 
of the compass. The prominent wind direction is as usual 
from the northwest -during the April to September upwelling 
period and shifts to the south in the first and last quarters , 
but this year there was a significantly high frequency of 
northwesterly winds during the last months of the year. The 
computer generated upwelling index (Figure 10) shows a pri- 
mary upwelling maximum in May, a second peak in about mid 
September and another in mid November. The possibility of 
a secondary upwelling peak is suggested by temperature obser- 
vations by Traganza et al . (1976). Figure 11 shows a rise 

in the isotherms peaking 30 to 45 days after the indicated 
wind initiated upwelling maximums in May and September. Al- 
though the water column does not respond instantaneously to 
surface wind stress, the delay noted here is excessive 
(Barham, 1957). The delay may be attributable to the assump- 
tion that the winds measured at Monterey Airport are actually 
representative of the winds over the bay, when in fact they 
may not be. The delay may also be caused by the lack of suf- 
ficient data to accurately depict the seasonal trends of 
temperature in the water column. 
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SIGNIFICANT WIND FIELD 
frequency vs direction 
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Figure 9. Frequency distribution of wind field direction at Monterey 
Peninsula Airport, 1974. Wind direction was corrected for 
topographic effects. 
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Figure 10. Upwelling index for Monterey from corrected airport wind 
observations, 1974. 
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Figure 11. Seasonal temperature variation off Monterey Bay. Values are 
the mean of four to nine stations (Traganza et al . , 1976). 



Traganza et al . (1976) observed a sharp rise in phosphate 

(Figure 12) with the major upwelling in May and a slight in- 
crease in phosphate and salinity (Figure 13) in September. 

No salinity data were available for the first half of the 
year. There was no clear-cut correlation of either phosphate 
or salinity to the upwelling index during November. Nutrient 
data from four to nine stations taken during seven cruises 
in 1974 were averaged (see Figure 12) and show a rapid rise 
in October which is probably not related to upwelling. From 
the five year study of Monterey Bay by Bolin (1964) it is 
known that the nutrient concentration is characteristically 
low over the depth of the euphotic zone (0-200 meters) for 
two or three months at the end of the year. The restoring 
of the nutrient level in all but the upper 20 to 30 meters 
to half of its May upwelling value while the surface tempera- 
ture reached a_ peak, may have signalled the beginning of the 
Davidson current period (see Bolin and Abbott, 1962). The 
Davidson current brings a southerly winter oceanic water 
mass into the Monterey region. This "Davidson water" is 
characterized by lower surface temperatures and surface salin- 
ities (due to high amounts of rain) . The deeper water of 
the euphotic zone has higher salinities and nutrient concen- 
trations than exist at Monterey. The mechanism which estab- 
lishes these characteristics may be winter storm mixing oc- 
curring south of Monterey or upwelling and mixing from below 
brought about by divergent (cyclonic) eddies formed in the 
current stream as it moves northward along the coast 
(Traganza et al., 1976). 
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Figure 12. Seasonal phosphate variation off Monterey Bay. Values 

are the mean of four to nine stations (Traganza et al., 1976). 



METERS 




Figure 13. Seasonal salinity variation off Monterey 
Bay. Values are the mean of four to nine 
stations (Traganza et al. , 1976). 
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Traganza's (1976) data were obtained at stations in an 
X pattern referenced to a drogue (Figure 14) and at inter- 
vals of about two miles, two and one-half hours apart. The 
data from four to nine stations were averaged to obtain the 
values shown in Figures 11, 12 and 13. 

Barham (1957) argues that the vertical circulation of 
Monterey Bay during the upwelling period is characterized 
by a region of surface divergence coincident with the head 
of the Monterey submarine canyon. The effect is due to 
"channeling" of upweiled replacement water along the canyon 
axis with a plume appearing at the surface near shore (see 
Figure 15). Evidence in support of this circulation concept 
is given by a distinct gradient of phytoplankton concentra- 
tion in the surface waters outlying the canyon head (Barham, 
1957). Barham believes that the gradient is due to the 
rapid horizontal advection of phytoplankton innoculum repre- 
senting a potential bloom away from the canyon head before 
the population shows any significant growth. During periods 
of reduced upwelling as indicated by the trends of the iso- 
therms , phytoplankton counts over the canyon rose to higher 
values , indicating that surface divergence or horizontal ad- 
vection had diminished. This possibility should not alter 
the timing of the model which is related to the general up- 
welling, but it does suggest serious spatial sampling prob- 
lems . 



53 




Figure 14- . Chart of the Monterey area showing station 
location . 
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Figure 15. Possible "fountainhead” effect over 
canyon head. 
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B. SENSITIVITY ANALYSIS 



The model was tested under various conditions of phyto- 
plankton sinking rates , predation pressure on herbivorous 
zooplankton, and advection of phytoplankton to determine the 
effect on the simulation of the seasonal cycle phosphate, 
phytoplankton and herbivorous zooplankton. 

Losses of phytoplankton by sinking and therefore food 
limitation on zooplankton was examined by setting the preda- 
tion terms on zooplankton to zero and varying the rate of 
sinking of algal cells, i.e. the availability of food. Sink- 
ing rates of zero, three, and six meters per day were used 
(Parsons and Takahashi , 1973 and Riley, 1965). The effects 
are shown in Figures 16, 17, and 18. Under conditions of 
zero sinking, a rapid rise in zooplankton "standing crop" 
in late summer reduces the phytoplankton biomass quickly to 
a minimum on day 215. The large decrease in phytoplankton 
allows the nutrient concentration to remain at a relatively 
high level (Figure 16) since the uptake of nutrients by 
phytoplankton is decreased while ongoing zooplankton excre- 
tion and mixing by upwelling add to the nutrient concentra- 
tion. The rate of growth of the herbivorous zooplankton 
population or "standing crop" is shown to decrease with pro- 
gressively higher phytoplankton sinking rates (Figure 18) 
as does the magnitude of the maximum zooplankton biomass 
reached during the year. The maximum also occurs later with 
increased sinking rate. 
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Figure 16. Simulated seasonal phosphate concentration with various phyto- 
plankton sinking rates; no predation. 
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Figure 17. Simulated seasonal phytoplankton concentration with various 
phytoplankton sinking rates; no. predation. 
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Figure 18 . Simulated seasonal zooplankton concentration with various 
phytoplankton sinking rates; no predation. 



The sinking rate was varied in a second test but preda- 
tor pressure was exerted on the zooplankton by employing 
the step function predator-prey relationship detailed earli- 
er. Additional predators were introduced in the simulation 

9 

at a zooplankton biomass level of 1.0 g C/m and were 

2 

removed when zooplankton biomass fell to 0.2 g C/m . Thresh- 
old values were determined empirically. The zooplankton bio- 

• • 2 2 

mass maximum is reduced from 1-10 g C/m to 0. 1-1.0 g C/m 

by the addition of predation pressure (Figure 21). The 
phosphate curves (Figure 19) show a marked response to phyto- 
plankton growth as evidenced by a lower nutrient minimum 
following the peak of the summer bloom. A considerable 
change in the zooplankton response, i.e. lower and earlier, 
is evidently due to reduction of food sources as phytoplank- 
ton is allowed to sink out of the mixed layer. The initial 
effect of increasing the sinking .rate from zero to three 
m/day on zooplankton is a shift in the occurrence of the peak 
to approximately 70 days later. Further increase of the 
sinking rate to six m/day results in a twofold decrease in 
the carbon biomass of zooplankton with an additional 20 day 
delay of the maximum. It is apparent that the rate of growth 
of zooplankton is slowed down and the maximum biomass occurs 
later because of food limitation but that predator pressure 
is responsible for limiting the magnitude of the zooplankton 
biomass. The zooplankton peaks also occur earlier when pre- 
dation is applied and also occur before the phytoplankton 
peaks (Figure 20). The simulation, therefore requires a 
positive sinking rate. 
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Figure 19. Simulated seasonal phosphate concentration with various phyto- 
plankton sinking rates; with predation. 
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Figure 20. Simulated seasonal phytoplankton concentration with various 
phytoplankton sinking rates; with predation. 
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Figure 21. Simulated seasonal zooplankton concentration with various 
phytoplankton sinking rates; with predation. 



The advection rate of phytoplankton was next varied 
while the sinking rate was held constant at 1.0 m/day 
(Bannister, 1974) and predator pressure was imposed in the 
same manner (and with the same thresholds) as in the second 
test. The results are shown in Figures 22, 23, and 24. A 
slower rate of growth of zooplankton is produced again when 
the advection coefficient, K, is increased from 0.0 to 0.3 
m 1 . Increasing the value of K has the direct effect of in- 
creasing the rate of phytoplankton advection. The single 
maximum of zooplankton simulated under conditions of high 
advection where K = 0 . 3 m ^ is due apparently to the sensi- 
tivity of the model's zooplankton growth equation to the 
mixed layer temperature and which permits a rapid growth 
of zooplankton during the mixed layer temperature minimum. 

The single peak of zooplankton (Figure 24) coincided with a 
relative temperature minimum about day 260. In Pearson's 
simulation a ten percent decrease in mixed layer temperature 
had a marked effect on zooplankton growth. This effect 
carried over to the partially modified simulation model. 

The effect of rapid zooplankton growth in a low temperature 
mixed layer, despite contradicting /actors of growth such 
as low food availability, is due to the sensitivity of the 
respiration term in the zooplankton growth equation to tem- 
perature. This condition is a peculiarity of the model and 
it is doubtful whether the real ecosystem behaves in the 
same manner. The zooplankton peak precedes the phytoplankton 
peak when advection is zero indicating the model needs an 
advection term. 
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Figure 22. Simulated seasonal phosphate concentration with various phytoplankton 
advection rates (with predation and algal sinking rate = 1.0 m/day). 
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1 igure 23. Simulated seasonal phytoplankton concentration with various phytoplankton 
advection rates (with predation and algal sinking rate = 1.0 m/day). 
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Figure 24. Simulated seasonal zooplankton concentration with various phytoplankton 
advection rates (with predation and algal sinking rate = 1.0 m/day). 



The predation pressure terms in the zooplankton equation 
were varied (Figures 25, 26, and 27) while keeping the sink- 
ing rate at 1.0 m/day and the advection coefficient at K = 

0.1 m ^ . The conditions of the test were: 

1. predation pressure set to zero; 

2. predation pressure defined by the linear function of 
zooplankton biomass: predator pressure = L times 
zooplankton biomass, where L = 0.01, and the units 
are g C/m^-day = (day -1 )(g C/m^ ) ; 

3. the linear function in 2, but L = 0.02; 

4. the step function discussed under the methods section 
of this thesis. 

In all cases, the resulting zooplankton growth appears food- 
limited until about day 120 when phytoplankton growth begins 
to show a rate change. Complete removal of predator pressure 
allows the zooplankton to increase rapidly until food limita- 
tion (be depletion of phytoplankton) again occurs on day 275. 

When the predator pressure increases linearly (with L = 0.01), 

2 

a single zooplankton peak of 3.05 g C/m occurs about day 
300. This is not consistent with observed conditions 
(Traganza, 1976). Doubling the predation pressure rate to 
L = 0.02 day ^ severely restricts zooplankton growth but pro- 
duces a curve with a hint of temporal conformity to the 
observed zooplankton maximum and minimum biomass levels. 

When the predator function is triggered at thresholds of 
zooplankton biomass, the zooplankton exhibit a double peak 
which is suggested by Traganza 's (1976) data. 
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0 100 200 300 

Figure 25. Simulated zooplankton response to various predation conditions (algal 
sinking rate = 1.0 m/day and advection coefficient = 0.1 m”-*-). 



C. COMPARISON OF SIMULATION RESULTS WITH OBSERVED DATA 



The computer simulation was run with the new upwelling 
and radiation indices and the step function predation pres- 
sure term detailed in the preceding section. The sinking 
rate was set at 1.6 m/day and the advection coefficient set 
to X = 0.1 m - '*". 

A comparison of the simulation results with data observed 
by Traganza et al. (1976) is shown in Figures 26, 27, and 28. 
The calculated nutrient .values lag the averaged mixed layer 
concentrations by 20 to 30 days during the summer months with 
a moderate amount of error in magnitude. The nutrient maxi- 
mum occurs on day 144 with a value of 2.08 pg-at P/1 in the 
simulation. A rapid decrease in phosphate levels coincident 
with the summer phytoplankton bloom displays the impact of 
nutrient uptake by phytoplankton in the model. 

The simulated phytoplankton peak occurs 40 days after 
the nutrient peak during the period of maximum incident radia- 
tion (see Appendix C). The modeled response of phytoplankton 
was relatively inflexible with regard to the radiation index, 
i.e. the phytoplankton peak was found to occur within a few 
days of the time of maximum incident radiation when the sink- 
ing rate was varied from about 1.0 m/day to 1.8 m/day. The 
relationship of the phytoplankton peak to the radiation peak 
is explained by examination of the phytoplankton growth equa- 
tions (shown in Appendix B). A slight increase in phytoplank- 
ton biomass was simulated during December (day 330) which 
could be due to the increased phosphate concentration at the 
end of the year. 
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Figure 26. Comparison of simulated seasonal phosphate with observed data. 
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Figure 27. Simulated seasonal phytoplankton concentration. 
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Figure 28 . Comparison of simulated seasonal zooplankton concentration with 
observed data. 



Twenty days after the phytoplankton maximum, the model 

simulated a zooplankton peak with a concentration of 0.97 
2 

g C/m . The zooplankton biomass growth corresponds well 

with the observed data until approximately the beginning of 

September. The measured values increase to a maximum of 
2 

1.85 g C/m in mid-December while the model simulated a 

2 

recovery peak of 0.47 g C/m and slightly earlier. 
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V. DISCUSSION 



The Monterey upwelling ecosystem simulation developed in 
this thesis has been shown to follow the seasonal trends of 
the observed phosphate, phytoplankton, and zooplankton data. 
There are errors in magnitude of the simulated response at 
various times of the year. The differences between observed 
and simulated zooplankton are most likely due to the fact 
that the model simulates herbivorous zooplankton while the 
observations include both herbivores and carnivores . Quanti- 
tative data defining the seasonal ratio of herbivorous to 
carnivorous zooplankton are needed to verify the simulation 
results. The encouraging aspect of the model is its ability 
to follow the seasonal trends and the fact that the timing 
of the response of phytoplankton biomass to the nutrient con- 
centration and the response of zooplankton to the phytoplank- 
ton cycle is biologically sound, i.e. a moderate delay is 
simulated between the peaks of successive trophic levels. 

Since all the factors likely to affect the three state 
variables (phosphate, phytoplankton and zooplankton) are not 
included in the model, some error in magnitude should appear. 
In recalling the philosophy expressed earlier, the objective 
of this thesis has been to create a time simulation of the 
dynamics of phosphate, phytoplankton and herbivorous zoo- 
plankton in a limited area over a long time. This objective 
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has been met by balancing the four model characteristics of 
generality, realism, precision and simplicity. 

From studies by Barham (1957) and Bolin (1964) it is 
noted that there may be characteristic seasonal patterns of 
plankton dynamics in response to characteristic patterns 
of the hydrography of the Monterey upwelling ecosystem. It 
is perhaps important to realize that the apparent limita- 
tions of the model as developed this far are relatively 
small when the use of the simulation to reproduce the 
characteristic patterns is considered per se . 

As more and better data with which to define the dynamics 
are obtained, the simulation model of the Monterey upwelling 
ecosystem can progress to a more refined level. 
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VI. SUGGESTIONS FOR FURTHER RESEARCH 



The following comments are made as a result of questions 
which arose during the modification of the Monterey ecosystem 
model. This list is provided to identify some areas which 
may warrant further investigation. It is by no means com- 
plete, but should serve as a guide. 

1. The most significant outcome of this research has 
been the recognition of the importance of accurate forcing 
functions. The presence of the submarine canyon in Monterey 
Bay and the likely effect on upwelling suggests further 
effort be directed to verifying the upwelling index. 

2. The Monterey Bay area is not homogeneous in physical, 
chemical and biological parameters as evidenced by the 
patchiness of ' biological samples and the observable spatial 
gradients of the physical and chemical properties. 

3. Additional work on the biological aspects of the 
model including verification and refinement of the predation 
terms and development of an appropriate kinetic expression 
is indicated. 

4. The simulation appears overly-sensitive to tempera- 
ture variations in the mixed layer. 

5. Various investigators have shown the importance of 
phytophagous fish in ecosystem models. A fish herbivore 
term should be investigated. 
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6. Some question as to the stability of the CSMP 360 
routine over long time simulation has been raised. This 
might be investigated and resolved by initializing the 
model part way through the year. 

7. The effect of varying the forcing functions in the 
revised model has not yet been studied. 

8. The advection expression might be expanded to in- 
clude the effects of all current systems influencing the 
Monterey region. 
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APPENDIX A: MONTEREY BAY UPWELLING PROGRAM 
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APPENDIX C: FORCING FUNCTIONS USED IN SIMULATION 
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